WEEK 9 SLIDES
Before we discuss HOW to build predictive policing systems, we need to ask: SHOULD we?
“A statistically ‘good’ model can still be socially harmful.”
Technical Questions:
Critical Questions:
We will learn the technical methods AND critically interrogate their use. Both skills are essential for ethical data science.
What vendors and police departments claim:
Sounds great, right?
. . .
But these claims rest on critical assumptions:
||–|–|| | 1st: Hotspots | Kernel Density | Past crime locations | KDE maps | | 2nd: Risk Terrain | Logistic Reg. | Crime + features | RTM software | | 3rd: ML | Random Forest, Neural Nets | Everything | PredPol, Palantir | | 4th: Person-Based | Network analysis | Social connections | Strategic Subject List |
Each generation claims to be more “objective” and “accurate”
But what if they’re all built on the same flawed foundation?
Traditional definition (data mining):
Extended definition (Richardson et al. 2019):
“Data derived from or influenced by corrupt, biased, and unlawful practices, including data that has been intentionally manipulated or ‘juked,’ as well as data that is distorted by individual and societal biases.”
1. Fabricated/Manipulated Data
2. Systematically Biased Data
3. Missing/Incomplete Data
4. Proxy Problems
From TV to Reality:
The Wire (2004): “If the crime rate doesn’t fall, you most certainly will”
Baltimore Reality (2008-2018):
Result: 55+ potential lawsuits, thousands of convictions questioned
Question: What happens when this data trains a predictive algorithm?
The Promise: Data-driven accountability, crime reduction
The Reality Revealed:
The Dual Strategy:
2013: Independent audit confirmed systematic data problems
The Confirmation Bias Loop:
False Positives (Predicted High Risk, No Crime Occurs)
Cost to individuals:
Cost to communities:
False Negatives (Predicted Low Risk, Crime Occurs)
Cost to individuals:
Cost to communities:
Who experiences each?
Typical Reforms Required:
What They DON’T Typically Address:
What it includes:
What community organizations noted is MISSING:
In 2019: CPD contracted with UChicago Crime Lab to “update” predictive approaches
Question: Using which data?
Before we build the model, we need to understand the methods
1 Setup & Data Preparation
2 Baseline Comparison
3 Feature Engineering
Using Abandoned Cars as “Disorder Indicator”:
4 Count Regression Models
5 Spatial Cross-Validation
6 Model Comparison
. . .
Does this actually predict crime? Or does it predict where certain communities call 311?
What assumptions are embedded in using “disorder” as a predictor?
Moran’s I (one value for entire study area)
Question answered: “Is there spatial clustering overall?”
Range:
Limitation: Doesn’t tell us WHERE clusters are
Example: “Chicago burglaries show positive spatial autocorrelation (I = 0.65)”
Local Moran’s I (one value per location)
Question answered: “Where are the clusters? Which specific locations?”
For each location:
Advantage: Maps showing WHERE patterns exist
Example: “South Side shows significant High-High clustering”
Before calculating spatial autocorrelation, we need to define “neighbors”
Spatial weights matrix (W):
Common approaches:
Queen Contiguity: Neighbors share edge OR vertex (corner)
[A][B][C]
[D][E][F]
[G][H][I]Cell E’s Queen neighbors: B, D, F, H (edges) + A, C, G, I (corners) = 8 neighbors
Rook Contiguity: Neighbors share only edges
Cell E’s Rook neighbors: B, D, F, H = 4 neighbors
Our fishnet grid uses Queen contiguity (most common for regular grids)
Raw weights: Binary (0 or 1)
Problem: Corner cells have fewer neighbors than interior cells
Solution: Row standardization (“W” style)
\[w^*_{ij} = \frac{w_{ij}}{\sum_j w_{ij}}\]
Effect: Each row sums to 1
Example:
Why? Makes cells comparable regardless of number of neighbors
\[I = \frac{n}{\sum_i \sum_j w_{ij}} \cdot \frac{\sum_i \sum_j w_{ij}(Y_i - \bar{Y})(Y_j - \bar{Y})}{\sum_i (Y_i - \bar{Y})^2}\]
Where:
Interpretation:
\[I_i = \frac{(Y_i - \bar{Y})}{\sum_j (Y_j - \bar{Y})^2 / n} \sum_j w_{ij}(Y_j - \bar{Y})\]
For each location \(i\):
Numerator: How different is location \(i\) from mean?
Denominator: Variance of all locations
Weight: How different are neighbors from mean?
Result: Local statistic for each location
Critical: Must test for statistical significance!
Visualizes relationship between:
Quadrant II | Quadrant I
(Low-High) | (High-High)
Outliers | HOTSPOTS
|
--|--
|
Quadrant III | Quadrant IV
(Low-Low) | (High-Low)
COLDSPOTS | OutliersInterpretation:
Problem: Not all clusters are “real” - could be random chance
Solution: Permutation test
Process:
Result: Only report clusters that are unlikely to occur by chance
|||-|–|-| | High-High | Positive | Above mean | Above mean | Hotspot | | Low-Low | Positive | Below mean | Below mean | Coldspot | | High-Low | Negative | Above mean | Below mean | Outlier (isolated high) | | Low-High | Negative | Below mean | Above mean | Outlier (isolated low) |
library(spdep)
# Step 1: Create spatial object
fishnet_sp <- as_Spatial(fishnet)
# Step 2: Define neighbors (Queen contiguity)
neighbors <- poly2nb(fishnet_sp, queen = TRUE)
# Step 3: Create spatial weights (row-standardized)
weights <- nb2listw(neighbors, style = "W", zero.policy = TRUE)
# Step 4: Calculate Local Moran's I
local_moran <- localmoran(
fishnet$abandoned_cars, # Variable of interest
weights, # Spatial weights
zero.policy = TRUE # Handle cells with no neighbors
)
# Step 5: Extract components
fishnet$local_I <- local_moran[, "Ii"] # Local I statistic
fishnet$p_value <- local_moran[, "Pr(z != E(Ii))"] # P-value
fishnet$z_score <- local_moran[, "Z.Ii"] # Z-score# Standardize the variable for quadrant classification
fishnet$standardized_value <- scale(fishnet$abandoned_cars)
# Calculate spatial lag (weighted mean of neighbors)
fishnet$spatial_lag <- lag.listw(weights, fishnet$abandoned_cars)
fishnet$standardized_lag <- scale(fishnet$spatial_lag)
# Identify High-High clusters
fishnet$hotspot <- 0 # Default: not a hotspot
# Criteria:
# 1. Value above mean (standardized > 0)
# 2. Neighbors above mean (spatial lag > 0)
# 3. Statistically significant (p < 0.05)
fishnet$hotspot[
fishnet$standardized_value > 0 &
fishnet$standardized_lag > 0 &
fishnet$p_value < 0.05
] <- 1
# Count hotspots
sum(fishnet$hotspot)
What you’ve identified:
Statistically significant clusters of high abandoned car reports
What this MIGHT mean:
Critical question:
Are these “high abandoned car” areas or “high reporting” areas?
For modeling:
We’re using this as a proxy for “disorder” but recognize it’s a flawed proxy
Our outcome: Burglary counts per grid cell (0, 1, 2, 3, …)
Linear regression assumes:
\[Y_i = \beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + \varepsilon_i\]
where \(\varepsilon_i \sim N(0, \sigma^2)\)
Problems for counts:
# Typical pattern for crime data
ggplot(fishnet, aes(x = countBurglaries)) +
geom_histogram(binwidth = 1, fill = "#440154FF", color = "white") +
labs(
title = "Distribution of Burglary Counts",
subtitle = "Most cells have 0-2 burglaries, few have many",
x = "Burglaries per Cell",
y = "Number of Cells"
) +
theme_minimal()Characteristics:
Appropriate for count data:
\[P(Y = k) = \frac{\lambda^k e^{-\lambda}}{k!}\]
Where:
Key property: Mean = Variance = \(\lambda\)
Examples:
Link \(\lambda\) to predictors via log link:
\[\log(\lambda_i) = \beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + ... + \beta_p X_{pi}\]
Equivalently:
\[\lambda_i = \exp(\beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + ... + \beta_p X_{pi})\]
Why log link?
On log scale:
\(\beta_1\) = change in log(expected count) per unit increase in \(X_1\)
On count scale (exponentiate):
\(\exp(\beta_1)\) = multiplicative effect on expected count
Examples:
|||-| | 0.14 | 1.15 | 15% increase per unit of X | | -0.22 | 0.80 | 20% decrease per unit of X | | 0.00 | 1.00 | No effect | | 0.69 | 2.00 | Doubling per unit of X |
# Fit Poisson model
model_poisson <- glm(
countBurglaries ~ Abandoned_Cars + Abandoned_Cars.nn + abandoned.isSig.dist,
data = fishnet,
family = poisson(link = "log")
)
# View results
summary(model_poisson)
# Exponentiate coefficients for interpretation
exp(coef(model_poisson))
# Example output:
# exp(coef)
# (Intercept) 0.234
# Abandoned_Cars 1.151
# Abandoned_Cars.nn 0.998
# abandoned.isSig.dist 0.999
# Interpretation:
# - Each additional abandoned car → 15.1% increase in expected burglaries
# - Each meter from nearest abandoned car → 0.2% decrease in expected burglariesPoisson assumption: Variance = Mean
Reality with crime data: Variance > Mean (often MUCH larger)
Why overdispersion occurs:
Check for overdispersion:
\[\text{Dispersion} = \frac{\text{Residual Deviance}}{\text{Degrees of Freedom}}\]
# Fit Poisson model
model_pois <- glm(
countBurglaries ~ Abandoned_Cars + Abandoned_Cars.nn + abandoned.isSig.dist,
data = fishnet,
family = poisson
)
# Calculate dispersion parameter
dispersion <- sum(residuals(model_pois, type = "pearson")^2) /
model_pois$df.residual
cat("Dispersion parameter:", round(dispersion, 3), "\n")
# Rule of thumb:
# < 1.5: Poisson OK
# 1.5 - 3: Mild overdispersion, NegBin recommended
# > 3: Serious overdispersion, NegBin essential
# Example output:
# Dispersion parameter: 4.523
# Warning: Serious overdispersion detected!Relaxes the variance = mean assumption
Adds dispersion parameter (\(\alpha\)):
\[\text{Var}(Y_i) = \mu_i + \alpha \mu_i^2\]
Where:
If \(\alpha = 0\): Reduces to Poisson
If \(\alpha > 0\): Allows extra variance (overdispersion)
Interpretation: Coefficients interpreted same way as Poisson!
library(MASS)
# Fit Negative Binomial model
model_nb <- glm.nb(
countBurglaries ~ Abandoned_Cars + Abandoned_Cars.nn + abandoned.isSig.dist,
data = fishnet
)
# View results
summary(model_nb)
# Compare to Poisson
AIC(model_pois) # e.g., 8234.5
AIC(model_nb) # e.g., 6721.3
# Lower AIC = better fit
# If NegBin AIC much lower → use NegBin
# Extract dispersion parameter (theta)
model_nb$theta # e.g., 2.47
# Interpretation: Significant overdispersion confirmed|–||-| | Variance assumption | Var = Mean | Var = μ + αμ² | | Overdispersion | Cannot handle | Accommodates | | Standard errors | Underestimated if overdispersed | Correctly estimated | | When to use | Count data, no overdispersion | Count data with overdispersion | | Crime data | Rarely appropriate | Usually better |
For today’s lab: We’ll fit both, compare, and use the better model
Unlike OLS, we don’t use residual plots the same way
Key diagnostics:
# Deviance residuals
plot(model_nb, which = 1) # Residuals vs. Fitted
# Identify outliers
outliers <- which(abs(residuals(model_nb, type = "deviance")) > 3)
# Influence
influence <- cooks.distance(model_nb)
influential <- which(influence > 4/length(influence))Basic model:
\[\log(\lambda_i) = \beta_0 + \beta_1 (\text{Abandoned Cars}_i)\]
Problem: Ignores spatial context
Enhanced model with spatial features:
\[\log(\lambda_i) = \beta_0 + \beta_1 (\text{Abandoned Cars}_i) + \beta_2 (\text{Distance to Nearest}_i) + \beta_3 (\text{Distance to Hotspot}_i)\]
Benefits:
Concept: Proximity matters for risk
Implementation: k-Nearest Neighbors (k=1)
For each grid cell:
In R:
library(FNN)
# Calculate distance to nearest abandoned car
nn_dist <- get.knnx(
data = st_coordinates(abandoned_cars), # "To" locations
query = st_coordinates(st_centroid(fishnet)), # "From" locations
k = 1 # Nearest 1
)
# Extract distances
fishnet$abandoned_car_nn <- nn_dist$nn.dist[, 1]Two-step process:
Step 1: Identify hotspots (Local Moran’s I High-High clusters)
Step 2: Calculate distance from each cell to nearest hotspot
# Step 1: Identify hotspots (we did this earlier)
hotspot_cells <- filter(fishnet, hotspot == 1)
# Step 2: Calculate distances
hotspot_dist <- get.knnx(
data = st_coordinates(st_centroid(hotspot_cells)),
query = st_coordinates(st_centroid(fishnet)),
k = 1
)
fishnet$hotspot_nn <- hotspot_dist$nn.dist[, 1]Why useful? Captures multi-scale spatial dependence:
# Create comparison maps
p1 <- ggplot(fishnet) +
geom_sf(aes(fill = abandoned_car_nn), color = NA) +
scale_fill_viridis_c(name = "Distance (m)", option = "plasma") +
labs(title = "Distance to Nearest Abandoned Car") +
theme_void()
p2 <- ggplot(fishnet) +
geom_sf(aes(fill = hotspot_nn), color = NA) +
scale_fill_viridis_c(name = "Distance (m)", option = "magma") +
labs(title = "Distance to Nearest Hotspot") +
theme_void()
grid.arrange(p1, p2, ncol = 2)Look for: Concentric patterns around features/hotspots
Problem: Crime is points, but we need areas for modeling
Options:
We use fishnet because:
library(sf)
# Step 1: Define cell size (in map units - meters for our projection)
cell_size <- 500 # 500m x 500m cells
# Step 2: Create grid over study area
fishnet <- st_make_grid(
chicago_boundary,
cellsize = cell_size,
square = TRUE,
what = "polygons" #you could change to hexagons if you wanted to be fancy.
) %>%
st_sf() %>%
mutate(uniqueID = row_number())
# Step 3: Clip to study area (remove cells outside boundary)
fishnet <- fishnet[chicago_boundary, ]
# Check results
nrow(fishnet) # Number of cells
st_area(fishnet[1, ]) # Area of one cell (should be 250,000 m²)
Common sizes:
Smaller cells:
Larger cells:
Choice affects results! No “correct” answer.
Process:
# Count burglaries per cell
burglary_counts <- st_join(burglaries, fishnet) %>%
st_drop_geometry() %>%
group_by(uniqueID) %>%
summarize(countBurglaries = n())
# Join back to fishnet
fishnet <- fishnet %>%
left_join(burglary_counts, by = "uniqueID") %>%
mutate(countBurglaries = replace_na(countBurglaries, 0))
# Summary
summary(fishnet$countBurglaries)
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 0 0 1 2.3 3 47Crime data typically has MANY zeros
Example distribution: - 40% of cells: 0 burglaries - 30% of cells: 1 burglary - 20% of cells: 2-3 burglaries - 10% of cells: 4+ burglaries
Implications:
Critical interpretation: Are zeros “true zeros” (no crime) or “missing data” (unreported)?
Standard k-fold CV:
Problem with spatial data:
Solution: Spatial cross-validation
Principle: Hold out entire spatial groups, not individual cells
Process:
Why better:
# Get unique districts
districts <- unique(fishnet$District)
# Initialize results
cv_results <- list()
# Loop through districts
for (dist in districts) {
# Split data
train_data <- fishnet %>% filter(District != dist)
test_data <- fishnet %>% filter(District == dist)
# Fit model on training data
model_cv <- glm.nb(
countBurglaries ~ Abandoned_Cars + Abandoned_Cars.nn + abandoned.isSig.dist,
data = train_data
)
# Predict on test data
test_data$prediction <- predict(model_cv, test_data, type = "response")
# Store results
cv_results[[dist]] <- test_data
}
# Combine all predictions
all_predictions <- bind_rows(cv_results)Common metrics for count models:
Mean Absolute Error (MAE): \[MAE = \frac{1}{n}\sum_i |y_i - \hat{y}_i|\]
Root Mean Squared Error (RMSE): \[RMSE = \sqrt{\frac{1}{n}\sum_i (y_i - \hat{y}_i)^2}\]
Mean Error (Bias): \[ME = \frac{1}{n}\sum_i (y_i - \hat{y}_i)\]
# Calculate metrics by district
cv_metrics <- all_predictions %>%
group_by(District) %>%
summarize(
MAE = mean(abs(countBurglaries - prediction)),
RMSE = sqrt(mean((countBurglaries - prediction)^2)),
ME = mean(countBurglaries - prediction)
)What do errors tell us?
MAE = 1.5 burglaries per cell:
RMSE = 2.3 burglaries per cell:
ME = -0.3 burglaries per cell:
Critical question: Are errors evenly distributed across neighborhoods?
# Map prediction errors
all_predictions <- all_predictions %>%
mutate(
error = countBurglaries - prediction,
abs_error = abs(error),
pct_error = (prediction - countBurglaries) / (countBurglaries + 1) * 100
)
# Visualize
ggplot(all_predictions) +
geom_sf(aes(fill = error), color = NA) +
scale_fill_gradient2(
low = "blue", mid = "white", high = "red",
midpoint = 0,
name = "Error"
) +
labs(title = "Prediction Errors",
subtitle = "Red = Over-prediction, Blue = Under-prediction") +
theme_void()Look for: Systematic patterns by neighborhood
KDE is the simplest spatial prediction method:
Benefits of comparison:
If our model doesn’t beat KDE, we’re wasting complexity
Conceptually:
Imagine placing a “hill” over each crime location:
Mathematical:
\[\hat{f}(s) = \frac{1}{n\tau^2}\sum_{i=1}^n K\left(\frac{d(s, s_i)}{\tau}\right)\]
Where:
Small bandwidth (e.g., 250m):
Large bandwidth (e.g., 2000m):
Common choices for urban crime:
Our approach:
library(spatstat)
# Step 1: Convert to point pattern (ppp) object
burglary_ppp <- as.ppp(
X = st_coordinates(burglaries),
W = as.owin(st_bbox(chicago_boundary))
)
# Step 2: Calculate KDE
kde_surface <- density.ppp(
burglary_ppp,
sigma = 1000, # Bandwidth in meters
edge = TRUE # Edge correction
)
# Step 3: Extract values to fishnet cells
fishnet$kde_risk <- raster::extract(
raster(kde_surface),
st_centroid(fishnet)
)
# Standardize to 0-1 scale for comparison
fishnet$kde_risk <- (fishnet$kde_risk - min(fishnet$kde_risk, na.rm=T)) /
(max(fishnet$kde_risk, na.rm=T) - min(fishnet$kde_risk, na.rm=T))Both KDE and our model produce continuous risk scores
For evaluation, convert to categories:
# Create quintiles (5 equal groups)
fishnet$model_risk_category <- cut(
fishnet$prediction,
breaks = quantile(fishnet$prediction, probs = seq(0, 1, 0.2)),
labels = c("1st (Lowest)", "2nd", "3rd", "4th", "5th (Highest)"),
include.lowest = TRUE
)
fishnet$kde_risk_category <- cut(
fishnet$kde_risk,
breaks = quantile(fishnet$kde_risk, probs = seq(0, 1, 0.2)),
labels = c("1st (Lowest)", "2nd", "3rd", "4th", "5th (Highest)"),
include.lowest = TRUE
)Each category contains 20% of cells
Critical: Test on data the model has NEVER seen
Process:
Question: Which method captures more crimes in high-risk areas?
# Spatial join: 2018 burglaries to fishnet with risk categories
results_2018 <- st_join(fishnet, burglaries_2018) %>%
group_by(model_risk_category) %>%
summarize(burglaries_2018 = n()) %>%
mutate(pct_of_total = 100 * burglaries_2018 / sum(burglaries_2018))Goal: Capture maximum crimes with minimum resources
Ideal result:
If we patrol the highest-risk 20% of cells (Category 5):
Interpretation:
Random patrol: Would capture ~20% (by chance)
Statistical Success:
✓ Model beats KDE
✓ Better than random
✓ Statistically significant difference
But ask:
A technically superior model can still be ethically problematic
Step 1: Data Preparation
Step 2: Exploratory Analysis
Step 3: Model Fitting
Step 4: Validation
Step 5: Comparison
Step 6: Critical Analysis
Instead of:
Consider:
Same technical methods, different normative goal
Examples:
Key differences: